In [1]:
import numpy as np
from numpy import random, fft, pi
import cft
from gnuplot import *
sigma = 1.0
N = 256
L = 22.0
B = cft.Box(2, N, L)
P = cft.Cutoff(B) * cft.Power_law(-1/2) * cft.Scale(B, 0.15)
S = cft.Scale(B, sigma)
peak = cft.make_peak(B, P, S)
def Hesse(B, P, S):
A = [cft.D(11), cft.D(22), cft.D(12)]
return [(a * S * cft.Potential())/(a * S * cft.Potential()).abs(B, P) for a in A]
hesse = Hesse(B, P, S)
pos = np.array([[0,0]]) + 11
#strength = np.r_[-3, -3] * 1 # , -1, -1] * 1.5
H = sum(([f * cft.Pos(x) for f in hesse] for x in pos), [])
h_clamp = [S * cft.Potential() * cft.Pos([11,11]) * cft.D(i) for i in [1,2]]
H += [h/h.abs(B, P) for h in h_clamp]
cgrf = cft.CGRF(B, P, H)
In [2]:
from IPython.display import Latex
def format_matrix(M):
s = "\\[\\left(\\begin{array}{" + "r"*M.shape[0] + "}\n "
s += " \\\\\n ".join([" & ".join(["{0: >6,.2f}".format(q) for q in row]) for row in np.array(M)])
s += "\n\\end{array}\\right)\\]\n"
return Latex(s)
format_matrix(cgrf.Q)
Out[2]:
In [3]:
G = np.array([-0.5, 1, 0, 0, 0]) * 3
#G = np.array([s * np.r_[1, 0, 0, alpha, alpha, 0] for s in strength]).flatten()
cgrf.generate_noise()
cgrf.set_coefficients(G)
A1 = cgrf.triplet()
A2 = cgrf.triplet(cft.Potential())
plot_six(B, *(A1 + A2))
print(cgrf.p_value())
In [4]:
from scipy.interpolate import RectBivariateSpline
x = np.mgrid[-B.L/2:B.L/2:B.N*1j]
Q = np.mgrid[-B.L/2:B.L/2:128j,-B.L/2:B.L/2:128j]
for i in range(3):
iPhi = RectBivariateSpline(x, x, A2[i])
sphi = iPhi.ev(Q[0].flat,Q[1].flat).reshape([128,128])
X = [(Q - np.gradient(D * sphi, B.L/128.)).transpose([1,2,0]).reshape([128**2, 2]) \
for D in [0.2, 0.8, 1.5]]
cmd = ("set term pngcairo size 900,300 font 'Bitstream Charter, 9'", "set multiplot layout 1,3",
"unset key", "set xrange [-10:10] ; set yrange [-10:10]") + \
tuple(reduce(lambda x,y: x+y, (["plot '-' w dots", x] for x in X))) + \
("unset multiplot",)
plot(*cmd)
In [4]:
In [5]:
from scipy.interpolate import RectBivariateSpline
from functools import reduce
x = np.mgrid[-B.L/2:B.L/2:B.N*1j]
Q = np.mgrid[-B.L/2:B.L/2:128j,-B.L/2:B.L/2:128j]
if False:
for i in range(3):
iPhi = RectBivariateSpline(x, x, A2[i])
sphi = iPhi.ev(Q[0].flat,Q[1].flat).reshape([128,128])
X = [(Q - np.gradient(D * sphi, B.L/128.)).transpose([1,2,0]).reshape([128**2, 2]) \
for D in [0.2, 0.5, 1.5]]
cmd = ("set term pngcairo size 720,240 font 'Bitstream Charter, 9'", "set multiplot layout 1,3",
"unset key", "set xrange [-10:10] ; set yrange [-10:10]") + \
tuple(reduce(lambda x,y: x+y, (["plot '-' w dots", x] for x in X))) + \
("unset multiplot",)
plot(*cmd)
Q = np.indices(B.shape) * B.res
for i in range(3):
vx = fft.ifftn(fft.fftn(A2[i]) * -1j * B.K[0]).real
vy = fft.ifftn(fft.fftn(A2[i]) * -1j * B.K[1]).real
X = (Q + 1.5 * np.array([vx, vy])).transpose([1,2,0]) + L
if (i == 1):
X += random.normal(0, 1e-5, X.shape)
T0 = X.reshape([B.N**2,2])
T1 = np.roll(X, -1, axis=0); T1[-1,:,0] += L ; T1 = T1.reshape([B.N**2,2])
T2 = np.roll(X, -1, axis=1); T2[:,-1,1] += L ; T2 = T2.reshape([B.N**2,2])
T3 = np.roll(np.roll(X, -1, axis=0), -1, axis=1)
T3[-1,:,0] += L ; T3[:,-1,1] += L ; T3 = T3.reshape([B.N**2, 2])
tt1 = np.c_[T0, T1, T2]
tt2 = np.c_[T3, T2, T1]
f = open('tmp{0}.dat'.format(i), 'wb')
f.write("{0}\n".format(2*len(tt1)).encode())
np.savetxt(f, tt1)
np.savetxt(f, tt2)
f.close()
In [6]:
import os
f = open('tmp.dat', 'wb')
for i in range(3):
os.system("../nerve2/nerve nerve -i tmp{0} -N 512 -L 22 --txt tmp{0}.dat".format(i))
np.savetxt(f, A1[i].T)
f.write("\n\n".encode())
f.close()
In [15]:
class Multiplot:
def __init__(self, w, h, n, m):
self.w = w ; self.h = h
self.n = n ; self.m = m
self.margin_left = 0.3
self.margin_bottom = 0.8
self.margin_top = 0.1
self.margin_right = 0.1
self.fw = w*n + self.margin_left + self.margin_right
self.fh = h*m + self.margin_top + self.margin_bottom
self.ml = self.margin_left / self.fw
self.mr = self.margin_right / self.fw
self.mt = self.margin_top / self.fh
self.mb = self.margin_bottom / self.fh
self.ew = (1 - (self.ml + self.mr)) / self.n
self.eh = (1 - (self.mt + self.mb)) / self.m
def set_term(self):
return "set term pdf size {0},{1} font 'Bitstream Charter, 10'".format(self.fw, self.fh)
def set_margins(self, i, j):
left = self.ml + self.ew * i
right = self.ml + self.ew * (i+1)
bottom = self.mb + self.eh * (self.m - j - 1)
top = self.mb + self.eh * (self.m - j)
return "set l{0} {1} ; set r{0} {2} ; set b{0} {3} ; set t{0} {4}".format(
"margin at screen", left, right, bottom, top)
def left(self, i):
return self.ml+self.ew*i
def bottom(self, j):
return self.mb + self.eh * (self.m - j - 1)
M = Multiplot(3, 3, 2, 3)
plot_pdf("filament-zeld",
M.set_term(),
"set xrange [-{0}:{0}] ; set yrange [-{0}:{0}]".format(L/2),
"set tics front",
"unset key", "unset colorbox", #"set size square",
"set multiplot",
"set cbrange [-2.5:2.5]",
"set xtics in format ''",
M.set_margins(0, 0),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 0),
M.set_margins(0, 1),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 1),
"set xtics format '%g'",
M.set_margins(0, 2),
"set colorbox horizontal user origin {0},{1}-{3} size {2}, {3}".format(M.left(0) + M.ew*0.05, M.bottom(2) - M.eh/10, M.ew*0.9, M.eh/20),
"plot 'tmp.dat' i {2} matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(N/L, L/2, 2),
"unset colorbox",
"set cbrange [0.02:5]", "set log cb",
"set ytics in format ''", "set xtics in format ''",
M.set_margins(1, 0),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 0),
M.set_margins(1, 1),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 1),
M.set_margins(1, 2),
"set xtics format '%g",
"set colorbox horizontal user origin {0},{1}-{3} size {2}, {3}".format(M.left(1) + M.ew*0.05, M.bottom(2) - M.eh/10, M.ew*0.9, M.eh/20),
"plot 'tmp{2}.nerve.conan' i 0 matrix u ($1/{0}-{1}):($2/{0}-{1}):3 t'' w image".format(512/L, L/2, 2),
"unset multiplot")
In [8]:
from scipy.interpolate import RectBivariateSpline
from thulsa import write_conan
f = open('mean.density.init.conan', 'wb')
write_conan(f, B, {"density": A1[1].T, "potential": A2[1].T})
f.close()
x = np.mgrid[-L/2:L/2:N*1j]
iPhi = RectBivariateSpline(x, x, A2[2])
In [9]:
os.system("../dmt2d/dmt2d lagrangian -I mean --truncate")
Out[9]:
In [10]:
gp = Gnuplot()
gp("set contour",
"set cntrparam levels incremental 0, 0.3, 3",
"set table", "set output 'gp_cntr.tab'",
"set view map", "unset surface",
# "splot 'mean.catastrophes.init.conan' i 0 matrix u ($1*%f - %f/2):($2*%f - %f/2):3" % (L/N, L, L/N, L),
"unset table")
gp.terminate() ; del gp
In [11]:
plot("set xrange [-3:3] ; set yrange [-2.5:2.5] ; set cbrange [-3:3]",
"set term pngcairo size 1000,700",
"set view map", "unset key",
"set size ratio 2.5/3",
"rcol(x) = 0.237 - 2.13*x + 26.92*x**2 - 65.5*x**3 + 63.5*x**4 - 22.36*x**5",
"gcol(x) = ((0.572 + 1.524*x - 1.811*x**2)/(1 - 0.291*x + 0.1574*x**2))**2",
"bcol(x) = 1/(1.579 - 4.03*x + 12.92*x**2 - 31.4*x**3 + 48.6*x**4 - 23.36*x**5)",
"set palette model RGB functions gcol(gray), rcol(gray), bcol(gray)",
"plot 'gp_cntr.tab' w l palette, \\",
" 'mean.a3.beta.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 4 lc rgb 'white', \\" % ((L/N, L)*2),
" 'mean.a3.alpha.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 4 lc rgb 'white', \\" % ((L/N, L)*2),
" 'mean.a3.beta.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 2 lc rgb '#332288', \\" % ((L/N, L)*2),
" 'mean.a3.alpha.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 2 lc rgb '#882255', \\" % ((L/N, L)*2),
" 'mean.points.init.conan' i 2 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 1.2 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 2 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 0.7 lc rgb '#88ccee', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 1 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 1.2 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 1 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 0.7 lc rgb '#cc6677', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 0 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 9 ps 1.8 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 0 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 9 ps 1.2 lc rgb 'black'" % (L/N, L, L/N, L))
In [12]:
from scipy.spatial import ConvexHull
pts = np.loadtxt("glass.txt")
pts = pts[np.where(abs(pts).max(axis=1) < L/2)]
v = iPhi.ev(pts[:,0], pts[:,1])
points = np.c_[pts[:,0], pts[:,1], (pts**2).sum(axis=1)/2 - 1.5 * v]
ch = ConvexHull(points)
data = ch.points[np.append(ch.simplices, ch.simplices[:,0:1], axis=1)]
f = open("tmp", 'w')
for x in data:
for y in x:
print(" ".join([str(l) for l in y]), file=f)
print("\n", file=f)
#gp = Gnuplot()
#gp("set term x11 font 'Bitstream Charter, 8'",
# "set xrange [-6:6] ; set yrange [-6:6]")
#gp("set hidden3d")
plot(#"set term pdf size 6,4 font 'Bitstream Charter, 10'", "set output 'tmp.pdf'",
"set term pngcairo size 900,600 font 'Bitstream Charter, 10'",
"set xrange [-5:5] ; set yrange [-3:3]", "unset key",
"plot 'tmp' w linespoints pt 7 ps 0.3 lw 0.3 lc rgb '#888888', \\",
" 'mean.a3.beta.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 4*2 lc rgb 'white', \\" % ((L/N, L)*2),
" 'mean.a3.alpha.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 4*2 lc rgb 'white', \\" % ((L/N, L)*2),
" 'mean.a3.beta.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 2*2 lc rgb '#332288', \\" % ((L/N, L)*2),
" 'mean.a3.alpha.init.conan' u ($1*%f - %f/2):($2*%f - %f/2):3 w l lw 2*2 lc rgb '#882255', \\" % ((L/N, L)*2),
" 'mean.points.init.conan' i 2 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 1.2*1 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 2 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 0.7*1 lc rgb '#88ccee', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 1 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 1.2*1 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 1 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 7 ps 0.7*1 lc rgb '#cc6677', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 0 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 9 ps 1.8*1 lc rgb 'white', \\" % (L/N, L, L/N, L),
" 'mean.points.init.conan' i 0 u ($1*%f - %f/2):($2*%f - %f/2) w p pt 9 ps 1.2*1 lc rgb 'black'" % (L/N, L, L/N, L))
#gp.terminate() ; del gp
In [13]:
def calc_normals(ch):
parab = np.zeros_like(ch.points)
parab[:,2] = (ch.points[:,:2]**2).sum(axis=1)
phi_pts = parab/2 - ch.points
def normalise(vec):
sign = vec[2]/abs(vec[2])
return vec/np.sqrt((vec**2).sum()) * sign
def flat(p):
return p[:2]
def normal(s):
if len(s) != 3:
return False
a, b, c = ch.points[s]
p, q, r = phi_pts[s]
n = np.cross(a - b, c - b)
pn = np.cross(p - q, r - q)
g = np.array([-pn[2]*pn[0], -pn[2]*pn[1], pn[0]**2 + pn[1]**2]) \
/ (pn[0]**2 + pn[1]**2)
A = np.cross(flat(a) - flat(b), flat(c) - flat(b))/2
p = (a + b + c)/3
return p, A, -n[:2]/n[2], g
return [normal(s) for s in ch.simplices if len(s) == 3 and ch.points[s[0]][2] < 16]
def vec2str(v):
return " ".join([str(i) for i in v])
V = calc_normals(ch) ; f = open('tmp', 'w')
for p, A, x, g in V:
print(vec2str(p), A, vec2str(x), vec2str(g), file=f)
#print(calc_normals(ch))
In [14]:
gp(#"set term pngcairo size 900,600", "set size ratio 2/3", "set xrange [-5:5] ; set yrange [-3:3]",
"set log cb", "set cbrange [0.01:1.5]", "set palette",
"plot 'tmp' u 1:2:(2*sqrt(abs($4))) t'' w p pt 7 ps variable,"
" 'tmp' u 1:2:($7/8):($8/8):(abs($4)) w vectors lc palette")
In [ ]:
gp = Gnuplot()
gp(#"set term pngcairo size 900,600", "set size 3/2",
#"set xrange [-5:5] ; set yrange [-3:3]",
"set palette ; set log cb ; set cbrange [0.01:2.0]",
"plot 'tmp' u 5:6:(5*sqrt(abs($4))) t'' w p pt 7 ps variable lc rgb '#aaffaa', \\",
" 'tmp' u 5:6:(5*sqrt(abs($4))) t'' w p pt 6 ps variable lc rgb '#44cc44', \\",
" 'tmp' u 5:6:($7/8):($8/8):(abs($4)) w vectors lw 1 lc palette")
# " '-' w p pt 7 ps 0.1 lc 3", X)
In [ ]: